%% generateDynamicEquation
% generate a 'dynamicEquation.m' file
%% Syntax
% generateDynamicEquation(Parameter)
%% Description
% Parameter: is a struct saving the model data
%
% generate a dynamicEquation.m file in the root directory


function generateDynamicEquation(Parameter)

% cheeck the running status parameters
if Parameter.Status.duration < 0 || Parameter.Status.acceleration < 0
    error('duration and acceleration must larger than 0 or equal 0');
end

if  Parameter.Status.vmax <= 0
    error('vmax must larger than 0');
end

if Parameter.Status.vmin > Parameter.Status.vmax
    error('vmin should smaller than vmax')
end

%% 

% check the exist of dynamic equation and create .txt
if isfile('dynamicEquation.m')
    delete dynamicEquation.m
end

dq = fopen('dynamicEquation.txt','w');


%%

% write comments, firstly

comments = [...
"%% dynamicEquation";...
"% saving the differential equation in this function";...
"%% Syntax";...
"% ddyn = dynamicEquation(tn, yn, dyn, Parameter)";...
"%% Description";...
"% tn: is the n-th time (s)";...
"% ";...
"% yn: is the displacement at the n-th time";...
"% ";...
"% dyn: is the first derivative at the n-th time";...
"% ";...
"% Parameter: is a struct saving all information about the model";...
"% ";...
"% ddyn: is the second derivative at the n-th time";...
" ";...
" "...
];

% function start
functionStart = [...
"function ddyn = dynamicEquation(tn, yn, dyn, Parameter)";...
" "...
];

fprintf(dq,'%s\n',comments);
fprintf(dq,'%s\n',functionStart);

%%

% calculate running status
Status   = Parameter.Status;
shaftNum = Parameter.Shaft.amount;
duration        = Status.duration;
ratio           = Status.ratio;
ratio           = [1; ratio]; % the first shaft is basic


vmax = zeros(1,shaftNum);
vmin = zeros(1,shaftNum);
acceleration = zeros(1,shaftNum);
for iShaft = 1:1:shaftNum
    vmax(iShaft)           = Status.vmax * ratio(iShaft);
    vmin(iShaft)           = Status.vmin * ratio(iShaft);
    acceleration(iShaft)   = Status.acceleration *ratio(iShaft);
end


%%

% write the code to calculate the omega 
% situation 1: a = 0
if (Status.acceleration == 0) && (Status.isDeceleration == false)
    
    calculateOmega = {...
     ' ';...
     '% calculate phase, speed and acceleration';...
    ['ddomega = zeros(1,',num2str(shaftNum),');'];...
    ['domega = [',num2str(vmax),'];'];...
     'omega = domega .* tn;';...
     ' '...
    };

elseif (Status.acceleration == 0) && (Status.isDeceleration == true)
    error('when acceleration equal 0, isDeceleration must be false')
end


% situation 2: a > 0 & no deceleration
if (Status.acceleration > 0) && (Status.isDeceleration == false)
    t1 = abs(vmax(1) / acceleration(1)); % the first key point (a>0 -> a=0)
    phase1 = 0.5 * acceleration .* t1^2; % the phase at t1
    
    calculateOmega = {...
     ' ';...
     '% calculate phase, speed and acceleration';...
    ['if tn <= ', num2str(t1)];...
    ['   ddomega = [', num2str(acceleration), '];' ];...
     '   domega  = ddomega .* tn ;';...
     '   omega   = 0.5 * ddomega * tn^2;';...
     'else';...
    ['   ddomega = zeros(1,', num2str(shaftNum),');'];...
    ['   domega  = [', num2str(vmax), '];'];...
    ['   omega   = [', num2str(phase1), '] + domega .* (tn-', num2str(t1), ');'];...
     'end';...
     ' '...
};

end


% situation 3: a > 0 & with deceleration
if (Status.acceleration > 0) && (Status.isDeceleration == true)
    t1 = abs(vmax(1) / acceleration(1)); % the first key point (a>0 -> a=0)
    phase1 = 0.5 * acceleration * t1^2; % the phase at t1
    t2 = t1 + duration; % the end of uniform motion (a=0 -> a<0)
    phase2 = phase1 + vmax * (t2 - t1);
    t3 = t2 + (abs(vmax(1))-abs(vmin(1)))/abs(acceleration(1)); % (a<0 -> a=0)
    phase3 = phase2 + vmax(1)*(t3-t2) - 0.5*acceleration*(t3-t2)^2;
    
calculateOmega = {...
 ' ';...
 '% calculate phase, speed and acceleration';...
['if tn <= ', num2str(t1)];...
['    ddomega = [', num2str(acceleration), '];'];...
['    domega  = [', num2str(acceleration),'] * tn;']
['    omega   = 0.5 * [', num2str(acceleration), '] * tn^2;'];...
['elseif tn <= ', num2str(t2)];...
['    ddomega = [', num2str(zeros(1,shaftNum)), '];'];...
['    domega  = [', num2str(vmax), '];'];...
['    omega   = [', num2str(phase1), '] + [', num2str(vmax), '] * (tn - ', num2str(t1), ' );'];...
['elseif tn <= ', num2str(t3)];...
['    ddomega = -[', num2str(acceleration), '];'];...
['    domega  = [', num2str(vmax), '] - [', num2str(acceleration), '] * (tn - ', num2str(t2), ' );'];...
['    omega   = [', num2str(phase2), '] + [', num2str(vmax), ']*(tn - ', num2str(t2), ' ) - 0.5*[', num2str(acceleration), ']*(tn - ', num2str(t2), ' )^2;'];...
'else';...
['    ddomega = [', num2str(zeros(1,shaftNum)), '];'];...
['    domega  = [', num2str(vmin), '];'];...
['    omega   = [', num2str(phase3), '] + [', num2str(vmin), ']*(tn - ', num2str(t3), ' );'];...
'end';...
' '
};
end


calculateOmega = cell2string(calculateOmega);
fprintf(dq,'%s\n', calculateOmega);

%%

% write load matrix part
loadMatrix1 = [...
" ";...
"% load matrix";...
"M = Parameter.Matrix.mass;";...
"G = Parameter.Matrix.gyroscopic;";...
"N = Parameter.Matrix.matrixN;";...
"Q = Parameter.Matrix.unblanceForce;"...
];
fprintf(dq,'%s\n', loadMatrix1);

% if there is loosing beaaring, the K and C would be created in bearingLoosingForce.m 
if ~Parameter.ComponentSwitch.hasLoosingBearing
    loadMatrix2 = [...
    "K = Parameter.Matrix.stiffness;";...
    "C = Parameter.Matrix.damping;";...
    ];
    fprintf(dq,'%s\n', loadMatrix2);
end


%%

% calculate the dof range of shaft
Node = Parameter.Mesh.Node;
dofInterval = Parameter.Mesh.dofInterval;

shaftDof = zeros(shaftNum,2);
for iShaft = 1:1:shaftNum
    IShaftNode = Node( [Node.onShaftNo] == iShaft & [Node.isBearing] == false );
    startNode = min([IShaftNode.name]);
    endNode = max([IShaftNode.name]);
    shaftDof(iShaft,:) = [dofInterval(startNode,1), dofInterval(endNode,2)];
end

%%

% write the part of processing G N Q
for iShaft = 1:1:shaftNum
processGN = {...
['G(', num2str(shaftDof(iShaft,1)), ':', num2str(shaftDof(iShaft,2)), ', ', num2str(shaftDof(iShaft,1)), ':', num2str(shaftDof(iShaft,2)), ')',...
' = domega(', num2str(iShaft), ')*G(', num2str(shaftDof(iShaft,1)), ':', num2str(shaftDof(iShaft,2)), ', ', num2str(shaftDof(iShaft,1)), ':', num2str(shaftDof(iShaft,2)), ');'];...
['N(', num2str(shaftDof(iShaft,1)), ':', num2str(shaftDof(iShaft,2)), ', ', num2str(shaftDof(iShaft,1)), ':', num2str(shaftDof(iShaft,2)), ')',...
' = ddomega(', num2str(iShaft), ')*N(', num2str(shaftDof(iShaft,1)), ':', num2str(shaftDof(iShaft,2)), ', ', num2str(shaftDof(iShaft,1)), ':', num2str(shaftDof(iShaft,2)), ');'];...
};
processGN = cell2string(processGN);
fprintf(dq,'%s\n', processGN);
end

%%

% calculate the dof range of disk
Disk = Parameter.Disk;
diskDof = dofInterval(Disk.positionOnShaftNode,:)';
diskDof = diskDof(1,:);
diskInShaftNo = Disk.inShaftNo';
diskMass = (pi * Disk.radius.^2 .* Disk.thickness) .* Disk.density;
diskMass = diskMass';
eccenticity = Disk.eccentricity';
coeffQ = (diskMass .* eccenticity);

%%

% write the part of unblance force
processQ = {...
 ' ';...
 ' ';...
 '% calculate unblance force';...
['diskInShaftNo = [', num2str(diskInShaftNo), '];'];...
['Q([', num2str(diskDof), '])   = [', num2str(coeffQ), '] .* ( ddomega(diskInShaftNo) .* sin(omega(diskInShaftNo)) + domega(diskInShaftNo).^2 .* cos(omega(diskInShaftNo)));'];...
['Q([', num2str(diskDof+1), ']) = [', num2str(coeffQ), '] .* (-ddomega(diskInShaftNo) .* cos(omega(diskInShaftNo)) + domega(diskInShaftNo).^2 .* sin(omega(diskInShaftNo)));'];...
 ' '...
};

processQ = cell2string(processQ);
fprintf(dq,'%s\n', processQ);

%%

% rub-impact force
if Parameter.ComponentSwitch.hasRubImpact
    generateRubImpactForce(Parameter.RubImpact, Parameter.Mesh)
    processRubImpactForce = [...
" ";...
"% calculate rub-impact force";...
"fRub = rubImpactForce(yn);";...
" ";...
];
    fprintf(dq,'%s\n', processRubImpactForce);
    rubForce = ' + fRub';
else
    rubForce = '';
end

%%

% loosing bearing

if Parameter.ComponentSwitch.hasLoosingBearing
    generateBearingLoosingForce(Parameter.LoosingBearing, Parameter.Mesh);
    proecessLoosingBearing = [...
" ";...
"% check the loosing bearing. If bearing loose, create the KLoose, CLoose.";...
"[K, C] = bearingLoosingForce(yn, Parameter.Matrix);";...        
" ";...
];
    fprintf(dq,'%s\n',  proecessLoosingBearing);
end

%%

% coupling misalignment force
if Parameter.ComponentSwitch.hasCouplingMisalignment
    generateMisalignmentForce(Parameter.CouplingMisalignment, Parameter.Mesh)
    processMisalignment = [...
" ";...
"% calculate couplingment misalignment force";...
"fMisalignment = misalignmentForce(omega, domega);";...
" ";...
];
    fprintf(dq,'%s\n', processMisalignment);
    misForce = ' + fMisalignment';
else
    misForce = '';
end

%%

% write total force
totalForce = {...
 ' ';...
 '% total force ';...
 ['F = Q', rubForce, misForce, ';'];...
 ' ';...
 };
totalForce = cell2string(totalForce);
fprintf(dq,'%s\n', totalForce);

%% 

% write dynamic equation
dynamicEquationStr = [...
" ";...
"% dynamic equation ";...
"ddyn = M \ ( F -  (K - N)*yn - (C - G)*dyn );";...
" ";...
];
fprintf(dq,'%s\n', dynamicEquationStr);

%%

% function end
functionEnd = [...
"end";...
" "...
];
fprintf(dq,'%s\n',functionEnd);

%%
% close .txt and transfer .txt -> .m
fclose(dq);
system('rename dynamicEquation.txt dynamicEquation.m');


end % end function